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ABSTRACT 

q Car is the only colliding-wind binary for which high-energy y rays are detected. Although 
the physical conditions in the shock region change on timescales of hours to days, the variabil¬ 
ity seen at GeV energies is weak and on significantly longer timescales. The y-ray spectrum 
exhibits two features that can be interpreted as emission from the shocks on either side of the 
contact discontinuity. Here we report on the first time-dependent modelling of the non-thermal 
emission in q Car. We find that emission from primary electrons is likely not responsible for 
the y-ray emission, but accelerated protons interacting with the dense wind material can ex¬ 
plain the observations. In our model, efficient acceleration is required at both shocks, with 
the primary side acting as a hadron calorimeter, whilst on the companion side acceleration 
is limited by the flow time out of the system, resulting in changing acceleration conditions. 
The system therefore represents a unique laboratory for the exploration of hadronic particle 
acceleration in non-relativistic shocks. 
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1 INTRODUCTION 

Particle acceleration up to very high energies in pulsar wind neb¬ 
ulae and supernova remnants is well established, with non-thermal 
emission seen from radio to TeV energies. The mechanism of diffu¬ 
sive shock acceleration (DSA; e.g. |Drury|1983[ > is usually invoked 
to explain the non-thermal emission in these systems, suggesting 
that all systems with strong shocks accelerate particles. The shocks 
present in some Galactic colliding wind binaries (CWBs) - binary 
systems with two massive stars and powerful winds - seem to sat¬ 
isfy all the criteria for particle acceleration and detectable high en¬ 
ergy emission: shock velocities of >1000km s~', available wind 
power of ~10 37 ergs>, and copious targets for the production of 
high-energy radiation: soft photon fields and high-density gas. In¬ 
deed non-thermal radio emission from some CWBs has been seen 
(e.g. |De Becker|2007) , and models predicting emission at y-ray en¬ 
ergies from such systems developed (e.g. Eichler & Usov|1993||Be-| 


|naglia & Romero|2003||Bednarek|2005| Reimer et al.|20067 ~bnly 

recently, however, has strong experimental evidence appeared for 
y-ray emission from such systems: in the unique case of q Car, us¬ 
ing AGILE ( |Tavani et al.|2009) and Fermi -LAT i Abdo et al. 2009jb 


,eyder| 


A hard X-ray tail is also seen from q Car j Viotti et al.|2004 
|et al.|2010| >. 

q Car is a binary system in a ~5.5-year orbit ^Damineli et al.| 
|2008j >. The masses and mass-loss rates of the stellar companions 
(Table |T|, together with the eccentricity of the CWB (e * 0.9), 
make this a very unusual system. The thermal X-ray emission asso¬ 
ciated with the wind collision region (WCR) and surrounding neb¬ 


ula has been extremely well studied (e.g. |Hamaguchi et al.|20 1 4| 
and references therein) and considerable theoretical work has gone 
into understanding this emission (e.g. |Pittard & Corcoran||2002| 
|Parkin et al.p0TT||Madura et al.|2013) . 

The LAT-detected emission above 200 MeV has been reported 
by|Abdo et al.| ( [20T0| >, |Famier et al. | ( |20 1 1 fr , and |Reitberger et ah] 
< |20 1 2) . The high-energy y-ray spectrum exhibits two distinct fea¬ 
tures: a low-energy component with a cutoff around 1 GeV, and a 
significantly harder component above ~ 10 GeV. Both components 
are found to be variable, but on longer timescales and less dra¬ 
matically than in X-rays. Upper limits from H.E.S.S. at energies 
above ~500 GeV are more restrictive than the extrapolated Fermi- 
LAT flux, implying a sudden drop in y-ray flux fAbramowski et al~] 
|2012[ >. The high gas densities present in the WCR of CWBs may 
result in the efficient production of 7r°-decay y rays (Farmer et al. 
[2011[ |Bednarek & Pabich|201 l| l, and lead to a calorimetric situa¬ 
tion where all energy injected into particle acceleration is radiated 
away on short timescales. In addition, the two shocks in this sys¬ 
tem have very different properties, and cannot be considered as a 
single system ( |Bednarek & Pabich|2011| >. The acceleration and in¬ 
teraction of non-thermal particles in q Car has been studied in pre¬ 
vious work (e.g. Tavani et al.|2009|[Abdo et al.|2010| [Farmer et al.| 


2011; Bed narek & Pabich|2011 1 , but the complex geometry, and 

all of the relevant phase-dependent timescales have so far not been 
fully considered. Here we present a 3D dynamical model combined 
with particle injection, propagation and interaction used to study 
the origin of y-ray emission from q Car. We consider the stellar 
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Table 1. Stellar parameters used throughout this work. 


Parameter 

Primary 

Companion 

Reference 

R,(Bo) 

100 

20 

1 

T,(K) 

25800 

30000 

2 

L,(10 6 Lo) 

4 

0.3 

2 

M(Mq yr -1 ) 

4.8 x 10~ 4 

1.4 x 10~ 5 

3 

Uoo (kms -1 ) 

500 

3000 

3 


References: (1) Hillier et ah] (2001) ; (2) |Davidson & Humphrey's] j 1997) ; 
(3) |Parkin et al.| ' 2009^ 


wind shocks of both stars in t] Car, and model the light curve and 
spectra from MeV to TeV energies. We first present the geometri¬ 
cal model and discuss the relevant timescales in the system before 
moving to the full model, results and discussion. 


2 MODEL 

2.1 Dynamical model 

To model the non-thermal emission from r\ Car we apply the dy¬ 
namical model introduced in Parkin & Pittard (20081. The orbit of 
the two stars is calculated in the centre-of-mass frame and the stel¬ 
lar winds are assumed to have reached terminal velocities before 
they interact. Orbital and stellar parameters are given in Tab. [I] 
Two regions are defined (Parkin & Pittard|2008) : 

The shock-cap is the region where the two stellar winds col¬ 
lide, and the flow accelerates outwards from the stagnation point 
along the contact discontinuity (CD). The ballistic point is defined 
where the flow speed reaches 85% of the primary terminal wind 
speed. We assume radial symmetry w.r.t. to the apex, but the orbital 
motion introduces a skew in the orientation of the shock-cap. The 
skew angle changes as a function of orbital phase i p, and can reach 
up to 38° but is well below 10° for 0.1 < ip < 0.9. The shape of 
the shock cap, tangential velocity of the flow and surface density is 
determined by momentum balance of the two winds ( jCanto et al.| 
[1996} >. For 100 p bins from </> = 0.0 to ip = 1.0 we generate a grid 
of points along the shock cap. This grid is formed by 36 azimuthal 
and 82 radial points. 

The ballistic flow is entered beyond the ballistic point, where 
the flow is assumed to be unaffected by the stars’ gravity, ram pres¬ 
sure of the winds or thermal pressure in the WCR. Beyond the bal¬ 
listic point, the wind material of the two stars is assumed to mix 
over a mixing length (see below) and to flow with a speed accord¬ 
ing to momentum balance of the two winds. 

A collapse of the WCR onto the companion star or dynami¬ 
cal instabilities in the WCR were proposed to explain the dramatic 
variability seen in thermal X-rays close to periastron (e.g.iCor coran| 
|2005| |Parkin et al.|2011) . The collapse is modelled by turning off 
particle acceleration in the inner 80% of the shock cap during the 
X-ray minimum (0.985 < ip < 1.025). 

Fig.[I|shows the overall geometry of the system, the motion of 
the two stars and the stagnation point of the two winds, the shock 
cap size, geometry at two phases, and the line-of-sight to the ob¬ 
server (projected onto the orbital plane). The asymmetric path of 
the stagnation point is due to the skew of the shock cap caused by 
the rapid motion of the companion around periastron. 



Figure 1. Geometry of the r\ Car system as a function of orbital phase over¬ 
laid on a transmissivity map at 250 GeV at ip = 0.04 in the plane of the two 
stars. The green solid, long-dashed and short-lines show the motions of the 
primary, companion and stagnation point around the centre-of-mass of the 
system. At ip = 0.04 and tp = 0.5 the positions of these elements in the sys¬ 
tem are represented respectively by closed and open symbols. The position 
of the shock cap is indicated by a blue line, the ballistic part is shown as 
azure segments that show the direction of movement of each of the compu¬ 
tational elements. The projected line of sight to the observer is shown with 
a thick black arrow (Madura et al.[20l'2j . 


2.2 Timescales 


Several different timescales in r] Car are relevant for particle ac¬ 
celeration and interaction. The acceleration time f acc of a particle 
in DSA is determined by the shock speed v s , the magnetic field 
strength B, and the diffusion in terms of the Bohm diffusion coeffi¬ 
cient K — tJ dCC Ks- face ~ 50 // acc (q! kins 1 ^^ a cc.GeV/) S. 

The magnetic field of the two stars in 77 Car is unknown. We 
adopt plausible surface magnetic field strengths of 100 G (see e.g. 
Walder et al. p0l2) , and spatial dependence following |Eichler ~&| 
Usov ( 1993), with a toroidal form for most of the orbit and varying 
between^. 1 < B < 10) G. 

The post-shock (PS) regions on both sides of the CD are very 
different in nature: along the entire orbit the primary shock is ex¬ 
pected to be radiative, whereas the companion shock is adiabatic. 
The slow dense primary wind collapses into a thin sheet of dense 
material, while the gas on the companion side flows out of the sys¬ 
tem without significant cooling. We calculate the thickness of the 
primary PS region for a radiative shock following jZhekov & Palla| 
( |2007) >. For the companion side, the thickness of the PS region can 
be calculated from mass conservation. For the wind parameters in 
Tab.[T| the typical thickness of the companion PS region (~1 AU) is 
orders of magnitude larger than the primary PS region (~10 ~ 4 AU). 
Note that densities in the WCR are typically 10 s cm~ 3 (10 7 cm~ 3 ) on 
the primary (companion) side at apastron, and an order of magni¬ 
tude higher close to periastron. 

Fig. 0 shows the key timescales as a function of particle en¬ 
ergy in the primary and companion PS regions for ip = 0.05. On 
both sides of the CD, electrons cool much faster than the typical 
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Figure 2. Key timescales in the primary (black) and companion (gray) 
shocked wind. The flow time of a gas packet to reach the shock cap edge 
is shown in long-dashed, the acceleration timescale in short-dashed, and 
the cooling timescales for electrons (t C ool,e> Coulomb, IC, and synchrotron 
losses) and protons (t coo i,p) in solid lines. r/ m athrmacc = 15 is assumed in 
both cases. 


time it takes gas to flow from the shock apex to the edge.Coulomb 
scattering dominates below lOOMeV, while inverse Compton (IC) 
scattering (E e ~ 10 GeV) and synchrotron cooling {E e > 100 GeV) 
dominate at higher electron energies. 

On the primary side the cooling time for protons is still shorter 
than the flow time. This suggests that for this shock an equilibrium 
is reached and that essentially no accelerated particles leave the 
shock cap. The cooling time of protons on the companion side, 
however, is longer than the flow time, which implies that protons 
are accelerated whilst moving along the shock cap and escape in to 
the ballistic flow. This more complex situation is treated separately 
in a semi-analytic model described below. 


2.3 Time-dependent model 

The dynamical model implies that the solid angles of the shock-cap 
as viewed from the primary and companion stars are 1.2 and 5.2 sr, 
respectively. In addition, only the wind kinetic power normal to the 
shock will be available for particle acceleration, limiting the avail¬ 
able power to 1.6(7. 6 ) x 10 36 ergs ' 1 on the primary (companion) 
side. We model the two components of emission detected by Fermi- 
LAT as arising from the two sides of the WCR as initially suggested 
by |Bednarek & Pabich| ( j20ll) . The low energy component, which 
is extremely luminous and has a cutoff around a few GeV, originates 
on the primary side, where the high PS density will limit particle 
acceleration and provide high emission efficiency through the in¬ 
teraction of all accelerated protons. The harder, fainter high-energy 
component originates on the companion side, where the lower den¬ 
sity will allow for acceleration limited only by the flow timescale, 
and result in lower luminosity due to most particles escaping with¬ 
out interacting. However, a certain level of mixing between the two 
layers of the WCR is needed to reach the detected emission levels, a 
phenomenon seen in hydrodynamical simulations of p Car ijParkm] 
|et al.| 201 1 | >. 

The faster cooling timescales of electrons (Fig. [2j, mean that 
much smaller values of 7 / acc are required than for protons, with 
super-Bohm acceleration needed for the high-energy component. 
Electron dominance of either component requires an electron to 
proton ratio of more than one (cf the commonly assumed 1 %) due 


to the additional electron emission produced down to MeV ener¬ 
gies. For these reasons we adopt the hadronic scenario in the fol¬ 
lowing. 

Below we give a general description of how time-dependent 
particle injection and y-ray production is treated in our model. We 
employ the radiative code used in Hinton & Aharonian (2007]). 


Primary side: Acceleration of particles in the PS region of the 
primary is in saturation and counterbalanced by losses. The max¬ 
imum energy a particle can reach depends on the exact location 
at which it enters the shock. Not only the magnetic field changes 
across the shock cap, but also the shock velocity, as the angle be¬ 
tween stellar wind and shock cap is changing. At the same time, 
radiation energy densities and gas density change. To calculate the 
emission from the primary side of the shock cap, we calculate the 
tangential velocity v, of the radiative layer for each annulus (|Canto| 


et al. 


1996). The shock velocity is given by v s = | ^jv 3 rim - ty. 


In¬ 
serting v s into t acc and solving for the energy where t acc = t coo i yields 
the maximum particle energy in each annulus. The power available 
for particle acceleration P availii = £ aV aii,i depends on the solid an¬ 
gle of the annulus, D, = 2;r(cos0,_i - cos0,), the primary mass- 
loss rate and wind speed: P ava jy = A constant fraction 

e p of the available wind power is assumed to go into accelerated 
protons P p j = r p P av aii,i • From this we calculate the CR proton in¬ 
jection spectrum and equilibrium y-ray spectrum in each annulus 
( [Ginzburg & Syrovatskii|1964||Zabalza et al.|20 lT} . Given the high 
densities in the primary wind and PS region, there is an additional 
contribution to the high-energy y-ray emission from charged pion 
decay and subsequent secondary electron emission. 

Companion side: As can be seen in Fig. [2] the lower PS 
density on the companion side results in accelerated protons los¬ 
ing only a small fraction of their energy to p-p collisions in the 
shock-cap. DSA therefore takes place under changing conditions 
as the relativistic particle population flows outwards in the shock- 
cap. To approximate this acceleration, we follow the CR population 
through each annulus outwards on the shock cap, and apply a semi- 
analytic acceleration scheme as follows. 

We follow the standard picture of DSA in non-relativistic 
shocks ( jBell|19781> : as particles enter the shock they are scattered 
from downstream to upstream by the turbulent wake and from up¬ 
stream to downstream by the Alfven waves generated by the ener¬ 
getic particles themselves attempting to escape upstream. After a 
time At, particles with initial energy of Eq will have a final energy 
In (E a ,/E 0 ) = t/f acc , where t acc = 3p. dCC K B vfr(r + 1 )/(r - 1), r is the 
shock compression ratio, and v s is the velocity of the shock. For 
each crossing there is a probability \u 2 lc, where un is the down¬ 
stream flow velocity, of the particles being advected downstream, 
resulting in a mean probability of remaining in the shock after a 
time At of lnfPAi) = -3 /(r - l)ln(E,/Eo). This process results 
in a downstream particle distribution with a power-law index of 
Pcr = —(r + 2)/(r— 1), which under strong shock conditions (r = 4) 
results in the well-known pc r = -2 index. 

We derive the relativistic particle distributions along the 
shock-cap by, at each time step At, injecting a fraction e a of the 
wind kinetic power in particles with energy E 0 = 1 GeV. The gain 
in energy due to acceleration over the time At is applied as above, 
and a fraction (1 - P&) of the particles at each given energy are 
lost downstream. On subsequent steps, only the particles remain¬ 
ing in the shock will continue to be accelerated, whereas particles 
lost downstream will be advected with the annulus until the ballistic 
flow is reached (note that very little energy is lost in the downstream 
region to p-p collisions, as can be seen in Fig. [2]). In addition to the 
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energy injected in particles at E 0 , additional power is required to 
provide the particle energy gain. For the properties of the shock 
on the companion side of r] Car, we have found that a total kinetic 
power fraction of e c ~ 15e° is used for particle acceleration. 

Towards the edge of the shock-cap, the reduced wind ram 
pressure owing to shock obliquity, and the increasing CR pressure 
from the particles being accelerated may cause a modification of 
the structure of the shock, leading to nonlinear effects in accelera¬ 
tion (see, e.g., [Maikov & O’C Drury|2001[ for a review). To account 
for this effect we adopt the semi-analytic approach of non-linear 
DSA by Berezhko & Ellison (1999), resulting in a hardening at the 
highest energies of up to peg. — -1.75. A deeper study of the ac¬ 
celeration process in CWB, including shock modification, will be 
discussed in a forthcoming paper ( jZabalza et al.|i n prep. I. 

Emission from the ballistic flow: Given the low gas density 
in the companion PS region, most of the accelerated protons will 
not interact, but leave the shock cap in the ballistic flow. Simu¬ 
lations indicate that the wind material is subject to numerous in¬ 
stabilities and that the radiative layer of gas mixes with the wind 
material of the two stars i jParkin et al. 12011) . We assume no mix¬ 
ing of the two stellar winds in the shock-cap region, but full mix¬ 
ing over a certain mixing length beyond the ballistic point. The 
two stellar winds and the radiative layer are assumed to form a 
region of mixed material at the edge of the shock cap at a dis¬ 
tance r 0 from the apex of the shock with thickness do and density 
Po = Wps.iPi + <7 ps ,2P2 + dwaliPwaliWo- In our model mixing oc¬ 
curs exponentially, over a characteristic scale equal to the shock 
cap radius. While moving away from the ballistic point, the density 
of the mixed material decreases as (r 0 /r) 2 . As the mixed material 
flows outwards the two stars orbit each other and a spiral structure 
of dense rings will form (cf. Fig. [TJ. The y-ray emission from p-p 
interactions in the ballistic flow is calculated in time steps much 
shorter than the interaction timescale, with the emission spectrum 
fixed to that found for p-p emission at the shock cap edge. We note 
that the diffusion timescale out of the flow region is always much 
longer than the flow timescale for the adopted value of rj, the mod¬ 
elled B-field, and the energy range considered. 

Pair-production absorption is significant for emission above 
~100GeV given the strong stellar radiation fields. We calculated 
the pair production opacity (e.g. |Dubus|2006} along the line of sight 
from each point on the shock cap and ballistic flow in each phase 
bin, considering the system geometry for the given phase and the 
density and direction of the primary and companion stellar radia¬ 
tion fields (with parameters as shown in Tab. [TJ. A transmission 
map for the orbital plane is shown in Fig.|T| showing how 250 GeV 
emission from behind the stars is completely suppressed. 


3 RESULTS 

Fig. 0 shows the y-ray spectrum for two phase bins including i) 
hadronic emission from accelerated protons on the primary and 
companion side, ii) hadronic emission from protons accelerated in 
the companion shock and interacting in the ballistic flow, and iii) 
emission from secondary electrons on the primary side. A fraction 
of €p,c = 20 % of the available wind power goes in to particle accel¬ 
eration at both primary and companion shocks in the curves shown 
and different values of r] are tested. 

In the y-ray lightcurve shown in Fig.[4]the low-energy compo¬ 
nent is dominated by the primary at almost all phases. Variability 
on the primary side is predicted only when the WCR collapses. 
For the rest of the orbit, the y-ray emission on the primary side 



Figure 4. y-ray lightcurve in the two energy bands used by |Reitberger et al*| 
|20I2) . Red lines show the primary contribution, green the (absorbed) com¬ 
panion contribution, and black the total predicted emission. ?7 = 15 and 
e = 0.2 are used. 

is constant due to the calorimetric behaviour and constant injected 
power, y-ray emission from the companion side is variable over the 
orbit due to the changing densities in the ballistic flow. The high- 
energy lightcurve is dominated by emission from protons acceler¬ 
ated on the companion side that escape and interact in the ballistic 
flow. Some residual y-ray emission in both components of the y- 
ray spectrum are expected even during the collapse of the WCR as 
a result of protons interacting in the ballistic flow that have been 
launched at earlier phases at the ballistic point (cf. Fig.[T|. 


4 DISCUSSION 

The model described above provides a reasonable level of agree¬ 
ment with the observed y-ray light curve and SED of rj Car. Fur¬ 
thermore, the broad features of the expected emission emerge from 
simple arguments based on the system geometry, mass flow and en¬ 
ergetics. Consideration of the time-dependent 3-D geometry of the 
system is, however, critical for modelling of the light-curve around 
periastron and the impact of pair-production absorption. We have 
shown that inclusion of a realistic geometry is important for ener¬ 
getics and the relative contribution of the two shocks and that the 
level of emission around 1 GeV requires extremely high-efficiency 
y-ray production, consistent (uniquely) with hadron calorimetry. 
We consider the dominance of the SED by emission from accel¬ 
erated protons and nuclei to be robust. Dominance by electrons of 
either component is very difficult due to energy-loss timescales and 
would certainly require an electron to proton ratio of > 1 . 

For the companion wind shock to reach the energies observed 
requires either fairly efficient acceleration with r] dcc ~ 10 , or much 
higher magnetic fields than the adopted stellar surface field of 
100 G. To reach the required flux levels at ~ 10 GeV there must 
be some mixing of the shocked and cooled primary wind with the 
PS flow of the companion, resulting in 0(10%) of accelerated pro¬ 
tons interacting. That mixing occurs on a characteristic scale of 
the order of the shock cap size seems plausible (see e.g. |Parkin| 
|et al.|20iT| . The fact that acceleration in this shock proceeds un¬ 
der changing conditions as particles flow around the shock cap is 
intriguing, and analogous to the situation in very young supernova 
remnants. Combined with the possibility of non-linear effects and 
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Figure 3. Model spectra and y-ray SED for the periastron (left) and the apastron phase bin as used in |Reitberger et al!]j2012| >. Red long-dashed lines indicate 
the primary contribution to the emission, assuming // acc = 15. Solid (long-dashed) blue lines show the observed companion emission for // acc = 5 (t/ acc = 15), 
and short-dashed thin gray lines show the intrinsic emission. Gray thick lines show the total emission for // acc = 5 (solid) and // acc = 15 (short-dashed). TeV 
data are from jAbramowski et al.|j2012} . 


subsequent spectral hardening in this system, this fact implies that 
77 Car may prove to be a very valuable system for testing our ideas 
about particle acceleration. These acceleration considerations are 
discussed in a forthcoming publication ( [Zabalza et al.|in prep.j l. 

The strong suppression of emission above ~100GeV due to 
pair-production absorption is unavoidable, given the constrained 
geometry and stellar luminosities. However, the fact that in our sce¬ 
nario the p-p emission is widely distributed over the mixing region 
of the ballistic flow changes significantly the phase-dependence of 
the absorption relative to the simplest assumptions. |Reitberger et ah| 
(2012) introduced an external X-ray absorber to explain the two- 
component y-ray spectrum, but there is no observational evidence 
for such an absorber. The energy density required to cause a spec¬ 
tral feature in the SED is orders of magnitude larger than the X-ray 
emission of any component in 77 Car ( [Hamaguchi et al.|2014[ . 

We are encouraged that with plausible assumptions and for the 
same adopted parameters e and 77 for the two shocks, the model pro¬ 
vides reasonable agreement with the measured data. Our scenario 
can be tested with better measurements around periastron, perhaps 
possible with the lifetime Fermi -LAT data and with HESS-II and 
CTA. Agreement with the measured light-curves is not possible 
without invoking some kind of collapse around periastron, and is 
strongly motivated by the X-ray observations. The fact that in 77 Car 
the highest energy particles mostly escape from the system, and that 
this situation is likely the same for other CWBs (with lower density 
winds) implies a contribution of CWBs to the Galactic CRs up to 
TeV energies. However, the fact that other known CWB systems 
typically have much lower mass-loss rates also implies that such 
systems are likely very difficult to detect in y rays. 
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